Cargamos el Dataset.

census<-read.csv("census2000_conLongitudLatitud.csv",sep=",",header = T,stringsAsFactors = T)
head(census)
##      ID      LocX     LocY RegDens RegPop MedHHInc MeanHHSz
## 1 00601 -66.74947 18.18010      70 19,143   $9,888     3.24
## 2 00602 -67.18025 18.36329      83 42,042  $11,384     3.10
## 3 00603 -67.13422 18.44862      86 55,592  $10,748     2.84
## 4 00604 -67.13700 18.49899      83  3,844  $31,199     3.00
## 5 00606 -66.95881 18.18215      65  6,449   $9,243     3.20
## 6 00610 -67.13605 18.28832      79 28,005  $12,629     3.01
str(census)
## 'data.frame':    33178 obs. of  7 variables:
##  $ ID      : Factor w/ 33178 levels "00601","00602",..: 1 2 3 4 5 6 7 8 9 10 ...
##  $ LocX    : num  -66.7 -67.2 -67.1 -67.1 -67 ...
##  $ LocY    : num  18.2 18.4 18.4 18.5 18.2 ...
##  $ RegDens : int  70 83 86 83 65 79 81 78 81 69 ...
##  $ RegPop  : Factor w/ 14798 levels "0","1","1,000",..: 3927 10542 12063 8031 12480 6944 13553 1185 5869 13806 ...
##  $ MedHHInc: Factor w/ 17817 levels "$0","$10,000",..: 17629 180 23 5013 17599 330 382 11 191 420 ...
##  $ MeanHHSz: num  3.24 3.1 2.84 3 3.2 3.01 2.84 2.97 2.97 2.64 ...

Convirtamos las variables en númericas.

library(dplyr)
census$MeanHHSz<-as.numeric(census$MeanHHSz)
census$RegPop<-as.numeric(gsub(",","",census$RegPop))
census$RegDens<-as.numeric(census$RegDens)
census$MedHHInc<- substring(census$MedHHInc,2)
census$MedHHInc<-as.numeric(gsub(",","",census$MedHHInc))
census$ID<-as.numeric(census$ID)
census<- census %>%
  filter(!is.na(ID))
head(census)
##   ID      LocX     LocY RegDens RegPop MedHHInc MeanHHSz
## 1  1 -66.74947 18.18010      70  19143     9888     3.24
## 2  2 -67.18025 18.36329      83  42042    11384     3.10
## 3  3 -67.13422 18.44862      86  55592    10748     2.84
## 4  4 -67.13700 18.49899      83   3844    31199     3.00
## 5  5 -66.95881 18.18215      65   6449     9243     3.20
## 6  6 -67.13605 18.28832      79  28005    12629     3.01

Analicemos los datos.

library(funModeling)
library(GGally)
profiling_num(census)
##   variable        mean      std_dev variation_coef       p_01        p_05
## 1       ID 16589.50000  9577.807952      0.5773416  332.77000  1659.85000
## 2     LocX   -91.08434    15.070689     -0.1654586 -132.85954  -121.19788
## 3     LocY    38.83039     5.359397      0.1380207   26.02384    29.95419
## 4  RegDens    50.50002    28.865519      0.5715943    2.00000     6.00000
## 5   RegPop  8596.97739 12978.758221      1.5096885    0.00000    62.00000
## 6 MedHHInc 38248.09386 17469.135891      0.4567322    0.00000 15625.00000
## 7 MeanHHSz     2.50071     0.595747      0.2382311    0.00000     1.79000
##          p_25        p_50        p_75        p_95        p_99      skewness
## 1  8295.25000 16589.50000 24883.75000 31519.15000 32846.23000  0.000000e+00
## 2   -97.21948   -88.30875   -80.38066   -72.80126   -69.38863 -1.265908e+00
## 3    35.38395    39.46048    42.10560    46.53666    48.77167  3.222690e-02
## 4    26.00000    51.00000    75.00000    95.00000    99.00000 -1.615443e-06
## 5   656.00000  2515.00000 11167.50000 36687.30000 55985.30000  2.349338e+00
## 6 28903.75000 35762.00000 45229.25000 69015.20000 97089.46000  1.470537e+00
## 7     2.36000     2.55000     2.74000     3.18150     3.87000 -1.889714e+00
##    kurtosis          iqr                   range_98                range_80
## 1  1.800000 16588.500000         [332.77, 32846.23]       [3318.7, 29860.3]
## 2  5.365383    16.838820  [-132.85954, -69.3886288] [-115.6307, -74.799528]
## 3  5.611973     6.721648 [26.02384267, 48.77167471] [32.0565574, 44.761814]
## 4  1.799749    49.000000                    [2, 99]                [11, 90]
## 5  9.896483 10511.500000               [0, 55985.3]          [177, 27304.9]
## 6 10.208649 16325.500000              [0, 97089.46]      [22211.1, 58394.6]
## 7 12.638813     0.380000                  [0, 3.87]            [2.13, 2.97]
status(census)
##   variable q_zeros    p_zeros q_na       p_na q_inf p_inf    type unique
## 1       ID       0 0.00000000    0 0.00000000     0     0 numeric  33178
## 2     LocX       0 0.00000000    0 0.00000000     0     0 numeric  32946
## 3     LocY       0 0.00000000    0 0.00000000     0     0 numeric  33141
## 4  RegDens       0 0.00000000 1013 0.03053228     0     0 numeric    100
## 5   RegPop    1013 0.03053228    0 0.00000000     0     0 numeric  14798
## 6 MedHHInc    1082 0.03261197    0 0.00000000     0     0 numeric  17817
## 7 MeanHHSz    1081 0.03258183    0 0.00000000     0     0 numeric    402
plot_num(census)

describe(census)
## census 
## 
##  7  Variables      33178  Observations
## --------------------------------------------------------------------------------
## ID 
##        n  missing distinct     Info     Mean      Gmd      .05      .10 
##    33178        0    33178        1    16590    11060     1660     3319 
##      .25      .50      .75      .90      .95 
##     8295    16590    24884    29860    31519 
## 
## lowest :     1     2     3     4     5, highest: 33174 33175 33176 33177 33178
## --------------------------------------------------------------------------------
## LocX 
##        n  missing distinct     Info     Mean      Gmd      .05      .10 
##    33178        0    32946        1   -91.08    16.12  -121.20  -115.63 
##      .25      .50      .75      .90      .95 
##   -97.22   -88.31   -80.38   -74.80   -72.80 
## 
## lowest : -176.63680 -174.19630 -171.70090 -170.40870 -170.27200
## highest:  -65.66116  -65.65988  -65.62762  -65.45604  -65.29258
## --------------------------------------------------------------------------------
## LocY 
##        n  missing distinct     Info     Mean      Gmd      .05      .10 
##    33178        0    33141        1    38.83    5.814    29.95    32.06 
##      .25      .50      .75      .90      .95 
##    35.38    39.46    42.11    44.76    46.54 
## 
## lowest : 17.96223 17.96453 17.97011 17.98167 17.98414
## highest: 70.13346 70.21520 70.47766 70.64090 71.29953
## --------------------------------------------------------------------------------
## RegDens 
##        n  missing distinct     Info     Mean      Gmd      .05      .10 
##    32165     1013      100        1     50.5    33.33        6       11 
##      .25      .50      .75      .90      .95 
##       26       51       75       90       95 
## 
## lowest :   1   2   3   4   5, highest:  96  97  98  99 100
## --------------------------------------------------------------------------------
## RegPop 
##        n  missing distinct     Info     Mean      Gmd      .05      .10 
##    33178        0    14798        1     8597    11725       62      177 
##      .25      .50      .75      .90      .95 
##      656     2515    11168    27305    36687 
## 
## lowest :      0      1      3      4      5, highest: 106415 108189 113935 114301 144024
## --------------------------------------------------------------------------------
## MedHHInc 
##        n  missing distinct     Info     Mean      Gmd      .05      .10 
##    33178        0    17817        1    38248    17811    15625    22211 
##      .25      .50      .75      .90      .95 
##    28904    35762    45229    58395    69015 
## 
## lowest :      0   2499   3750   4318   4444, highest: 183686 185466 196298 200000 200001
## --------------------------------------------------------------------------------
## MeanHHSz 
##        n  missing distinct     Info     Mean      Gmd      .05      .10 
##    33178        0      402        1    2.501    0.529    1.790    2.130 
##      .25      .50      .75      .90      .95 
##    2.360    2.550    2.740    2.970    3.181 
## 
## lowest : 0.00 0.33 0.50 0.67 0.71, highest: 6.40 7.18 7.31 8.33 8.49
## --------------------------------------------------------------------------------

Podemos ver algunos valores característicos de los datos así como su distribución.

Alta presencia de NA,s en REgDens, presencia de 0 en RegPop y MeanHHSD, de hecho, RegDens tiene el mismo número de NA,s que que ceros REgPop, piede ser que coincidan. Veamos los NA,s gráficamente.

library(VIM)
aggr(census, col=c('navyblue','yellow'),
     numbers=TRUE, sortVars=TRUE,
     labels=names(census), cex.axis=.7,
     gap=3, ylab=c("Missing data","Pattern"))

## 
##  Variables sorted by number of missings: 
##  Variable      Count
##   RegDens 0.03053228
##        ID 0.00000000
##      LocX 0.00000000
##      LocY 0.00000000
##    RegPop 0.00000000
##  MedHHInc 0.00000000
##  MeanHHSz 0.00000000

Veamos si los NA,s de RegDens coinciden con los ceros de RegPop y de MeanHHsz.

census %>%
  filter(is.na(RegDens))%>%
  group_by(RegDens)%>%
  summarise(RegPop=unique(RegPop),MeanHHSz=unique(MeanHHSz))
## # A tibble: 1 x 3
##   RegDens RegPop MeanHHSz
##     <dbl>  <dbl>    <dbl>
## 1      NA      0        0

Todos los valores NA,s de RegDens coinciden con cers en las otras dos variables, por tanto eliminamos las filas con NA,s

census<-na.omit(census)

Discretizamos las variables.

library(nima)
census.scale<-census
census.scale$discreteRegDens <- discrete_by_quantile(census.scale$RegDens)/4
census.scale$discreteRegPop <- discrete_by_quantile(census.scale$RegPop)/4
census.scale$discreteMedHHInc <- discrete_by_quantile(census.scale$MedHHInc)/4
census.scale$discreteMeanHHSz <- discrete_by_quantile(census.scale$MeanHHSz)/4
summary(census.scale)
##        ID             LocX              LocY          RegDens     
##  Min.   :    1   Min.   :-176.64   Min.   :17.96   Min.   :  1.0  
##  1st Qu.: 8328   1st Qu.: -97.22   1st Qu.:35.40   1st Qu.: 26.0  
##  Median :16635   Median : -88.35   Median :39.47   Median : 51.0  
##  Mean   :16610   Mean   : -91.10   Mean   :38.85   Mean   : 50.5  
##  3rd Qu.:24880   3rd Qu.: -80.41   3rd Qu.:42.11   3rd Qu.: 75.0  
##  Max.   :33178   Max.   : -65.29   Max.   :71.30   Max.   :100.0  
##      RegPop          MedHHInc         MeanHHSz     discreteRegDens
##  Min.   :     1   Min.   :     0   Min.   :0.000   Min.   :0.25   
##  1st Qu.:   755   1st Qu.: 29632   1st Qu.:2.380   1st Qu.:0.25   
##  Median :  2720   Median : 36154   Median :2.560   Median :0.50   
##  Mean   :  8868   Mean   : 39453   Mean   :2.579   Mean   :0.62   
##  3rd Qu.: 11714   3rd Qu.: 45641   3rd Qu.:2.750   3rd Qu.:0.75   
##  Max.   :144024   Max.   :200001   Max.   :8.490   Max.   :1.00   
##  discreteRegPop   discreteMedHHInc discreteMeanHHSz
##  Min.   :0.2500   Min.   :0.250    Min.   :0.2500  
##  1st Qu.:0.2500   1st Qu.:0.250    1st Qu.:0.2500  
##  Median :0.5000   Median :0.500    Median :0.5000  
##  Mean   :0.6249   Mean   :0.625    Mean   :0.6191  
##  3rd Qu.:0.7500   3rd Qu.:0.750    3rd Qu.:0.7500  
##  Max.   :1.0000   Max.   :1.000    Max.   :1.0000
colnames(census.scale)[1]<-"ID"

Comprobemos la correlación entre las nuevas variables.

matrizCorrelacion<-cor(census.scale[,8:11],method=c("spearman"))
library(corrplot)
corrplot(matrizCorrelacion,method="number",type="upper")

No observamos una correlación excesivamente alta.

Para hacer clustering con el método jerárquico, primero extraeremos una muestra del 30%.

set.seed(2354)
indexes<-sample(1:nrow(census.scale),size=0.3*nrow(census.scale))
census.scale.muestra<-census.scale[indexes,]
dim(census.scale.muestra)
## [1] 9649   11

Hacemos la matriz de distancias y dibujamos un dendograma.

memory.limit(size=40000)
## [1] 40000
library(vegan)
library(cluster)
MatrizDistancias<-vegdist(census.scale.muestra[,8:11],method="euclidean")
clusterJerarquico<-hclust(MatrizDistancias,method="ward.D2")
plot(clusterJerarquico, labels = FALSE, main = "Dendrograma")
rect.hclust(clusterJerarquico, k=2, border="red") 
rect.hclust(clusterJerarquico, k=3, border="blue") 
rect.hclust(clusterJerarquico, k=4, border="green") 
rect.hclust(clusterJerarquico, k=5, border="yellow") 
rect.hclust(clusterJerarquico, k=6, border="purple") 
rect.hclust(clusterJerarquico, k=7, border="gray") 
rect.hclust(clusterJerarquico, k=8, border="black") 

A primera vista, a medida que va aumentando los clusters la ganancia disminuye, aún así, es difícil ver con exactitud cuantos clusters coger.

Veamos el método calinski y Harabasz. Este método escoge el número de clusters que maximizar la suma de las distancias entre clusters minimizando las distancias de los valores que pertenecen a cada cluster.

calinsky <- cascadeKM(census.scale.muestra[,8:11], inf.gr = 2, sup.gr = 10, iter = 100, criterion = "calinski")
calinsky$results
##          2 groups 3 groups 4 groups 5 groups 6 groups  7 groups  8 groups
## SSE      1937.370 1603.104 1333.785 1178.314 1042.926  944.5641  858.7526
## calinski 5298.448 4206.936 4019.761 3730.373 3621.717 3499.3864 3436.4828
##           9 groups 10 groups
## SSE       777.6642  714.2258
## calinski 3445.7626 3429.7318

Ambos métodos marcan que lo más óptimo es escoger dos clusters.

Veamos el gráfico de silueta, ponemos como máximo 20 clusters ya que no cogeremos más bajo ninguna circunstancia, los resultados deben ser interpretables.

asw <- numeric(20)
for(k in 2:(20 - 1)){
  sil <- silhouette(cutree(clusterJerarquico, k = k), MatrizDistancias)
  asw[k] <- summary(sil)$avg.width}
k.best <- which.max(asw)

plot(1: 20, asw, type="h", 
     main = "Silhouette-optimal number of clusters", 
     xlab = "k (number of groups)", ylab = "Average silhouette width")
axis(1, k.best, paste("optimum", k.best, sep = "\n"), col = "red", font = 2,
     col.axis = "red")
points(k.best, max(asw), pch = 16, col = "red", cex = 1.5)

Nos indica que debríamos coger dos clusters, según este método.

Apliquemos otro método.

library(GGally)
library(factoextra)
fviz_nbclust(census.scale.muestra, kmeans, method = "silhouette") +
  labs(subtitle = "Silhouette method")

El método Silhouette mide la calidad de una agrupación y determina cómo se encuentra cada punto dentro de su agrupación. Indica que los más óptimo es coger 3 clusters.

Se aplica el gráfico de Elbow al Dataset completo, aunque sea un método propio de k-means, nos puede ayudar a saber cuantos clusters coger, mi intención no es coger sólo 2, ya que pretendo ver un mayor grado de diferenciación entre los datos. Seguimos el mismo criterio, no cogeremos más de 20 clusters.

set.seed(13345)
n <- dim(census.scale)[1] 

p <- dim(census.scale[,8:11])[2]

SSW <- (n - 1) * sum(apply(census.scale[,8:11],2,var)) 

for (i in 2:20) SSW[i] <- 
  sum(kmeans(census.scale[,8:11],centers=i,nstart=3,iter.max=20)$withinss)

plot(1:20, SSW, type="b", xlab="Number of Clusters", ylab="Sum of squares within groups",pch=19, col="steelblue4")

Veamos el resultado de la gráfica numéricamente.

mejora<-numeric()
anterior<-1
x<-1
for (i in SSW){
  if (anterior==1){mejora[x]<-1}
  else {mejora[x]<-(anterior-i)/anterior}
  anterior <- i
  x <- x+1
}
round(mejora,3)
##  [1] 1.000 0.353 0.168 0.169 0.120 0.091 0.119 0.084 0.085 0.056 0.079 0.041
## [13] 0.032 0.041 0.079 0.012 0.036 0.061 0.000 0.014

Observando la evolución de la suma de cuadrado intragrupos observo que el porcentje de explicación de las variables empieza a aumentar en menos del 10% a partir del séptimo cluster. Tomaremos este valor, no quiero aumentar más el número de clusters para no perder demasiada interpretabilidad.

Calculemos los 7 centroides.

asignacionJerarquica<-cbind(census.scale.muestra[,c(1:3,8:11)],
                            cutree(clusterJerarquico, k = 7))
colnames(asignacionJerarquica)[8]<-"cluster"
centroidesJerarquico<-asignacionJerarquica%>%
  group_by(cluster)%>%
  summarise(total=n(),discreteRegDens=mean(discreteRegDens),discreteRegPop=mean(discreteRegPop),
            discreteHHInc=mean(discreteMedHHInc),discreteMeanHHSz=mean(discreteMeanHHSz))
## `summarise()` ungrouping output (override with `.groups` argument)
centroidesJerarquico
## # A tibble: 7 x 6
##   cluster total discreteRegDens discreteRegPop discreteHHInc discreteMeanHHSz
##     <int> <int>           <dbl>          <dbl>         <dbl>            <dbl>
## 1       1   951           0.357          0.488         0.358            0.866
## 2       2  2072           0.369          0.480         0.372            0.397
## 3       3  1772           0.906          0.926         0.674            0.402
## 4       4   616           0.846          0.317         0.492            0.435
## 5       5  1924           0.872          0.889         0.834            0.894
## 6       6  1286           0.402          0.460         0.760            0.557
## 7       7  1028           0.557          0.446         0.848            0.912

Observemos que porcentaje de nuestros datos representa cada cluster.

porcentaje<-cbind(centroidesJerarquico$cluster,porcentaje=centroidesJerarquico$total/
                    sum(centroidesJerarquico$total))
porcentaje
##        porcentaje
## [1,] 1 0.09855944
## [2,] 2 0.21473728
## [3,] 3 0.18364597
## [4,] 4 0.06384081
## [5,] 5 0.19939890
## [6,] 6 0.13327806
## [7,] 7 0.10653954

Introdujamos los centroides en el método k-means para poder ver así los clusters definitivos.

kmeans <- kmeans(census.scale[,8:11],centers=centroidesJerarquico[,3:6])
kmeans$centers
##   discreteRegDens discreteRegPop discreteMedHHInc discreteMeanHHSz
## 1       0.4011372      0.4770248        0.3735786        0.8629612
## 2       0.3454991      0.4532027        0.3695585        0.3803638
## 3       0.8938675      0.9282639        0.7053681        0.3915393
## 4       0.8745450      0.4082120        0.3612375        0.4419927
## 5       0.8752535      0.8920658        0.8341005        0.8911297
## 6       0.4534075      0.4268081        0.8224096        0.3916029
## 7       0.4867647      0.4719412        0.8460588        0.8731765

Estos son nuestros clusters definitivos, los introducimos en nuestro Dataset.

census<-cbind(census,kmeans$cluster)
head(census)
##   ID      LocX     LocY RegDens RegPop MedHHInc MeanHHSz kmeans$cluster
## 1  1 -66.74947 18.18010      70  19143     9888     3.24              5
## 2  2 -67.18025 18.36329      83  42042    11384     3.10              5
## 3  3 -67.13422 18.44862      86  55592    10748     2.84              5
## 4  4 -67.13700 18.49899      83   3844    31199     3.00              5
## 5  5 -66.95881 18.18215      65   6449     9243     3.20              1
## 6  6 -67.13605 18.28832      79  28005    12629     3.01              5
colnames(census)[8]<-"cluster"

Hagamos los gráficos de radar para ver los valores de cada cluster.

library(fmsb)



centroidesOptimizacionParaRadar<-rbind(
  rep(1,4) , 
  rep(0,4) , 
  apply(kmeans$centers , 2, mean),
  kmeans$centers)

colors_border=c( rgb(0.2,0.5,0.5,0.9), rgb(0.8,0.2,0.5,0.9) , rgb(0.7,0.5,0.1,0.9) )
colors_in=c( rgb(0.2,0.5,0.5,0.4), rgb(0.8,0.2,0.5,0.4) , rgb(0.7,0.5,0.1,0.4) )


  tamanyo<-centroidesJerarquico[1,2]
  
  radarchart( as.data.frame(centroidesOptimizacionParaRadar[c(1:3,4),])  , axistype=1 ,
              pcol=colors_border , pfcol=colors_in , plwd=1 , plty=1,
              cglcol="grey", cglty=1, axislabcol="grey", caxislabels=seq(0,1,5), cglwd=0.8,
              vlcex=0.8,
              title=paste0("Tamaño:",tamanyo))

Los habitantes del primer grupo forman parte de un area con familias muy grandes, la densisdad está por debajo de la media, el número de habitantes es ligeramente menor a la media, pueden ser ciudades grandes con baja densidad, los ingresos están por debajo de la media, llamaremos a este cluster “nucleos familiares grandes en ciudades poco densas”.

  tamanyo<-centroidesJerarquico[2,2]
  
  radarchart( as.data.frame(centroidesOptimizacionParaRadar[c(1:3,5),])  , axistype=1 ,
              pcol=colors_border , pfcol=colors_in , plwd=1 , plty=1,
              cglcol="grey", cglty=1, axislabcol="grey", caxislabels=seq(0,1,5), cglwd=0.8,
              vlcex=0.8,
              title=paste0("Tamaño:",tamanyo))

Este cluster cuenta con una densidad ligeramente por debajo de la media y pocas personas por región, pueden ser areas rurales, el tamaño familiar es bajo y los ingresos también, lo llamaremos “areas rurales empobrecidas”.Cabe destacar que este es el cluster con un mayor número de pertenecientes.

  tamanyo<-centroidesJerarquico[3,2]
  
  radarchart( as.data.frame(centroidesOptimizacionParaRadar[c(1:3,6),])  , axistype=1 ,
              pcol=colors_border , pfcol=colors_in , plwd=1 , plty=1,
              cglcol="grey", cglty=1, axislabcol="grey", caxislabels=seq(0,1,5), cglwd=0.8,
              vlcex=0.8,
              title=paste0("Tamaño:",tamanyo))

Este cluster cuenta con una densidad de población y número de habitantes muy alto, deben ser ciudades importantes, los ingresos están ligeramente por encima de la media mientras que el tamaño familiar es bajo, por lo tanto puden ser areas donde la gente se muda para trabajar, posiblemente muchas personas de este cluster sean jóvenes que se han independizado y que aún no han formado una familia, lo llamaremos “grandes ciudades juveniles”.

  tamanyo<-centroidesJerarquico[4,2]
  
  radarchart( as.data.frame(centroidesOptimizacionParaRadar[c(1:3,7),])  , axistype=1 ,
              pcol=colors_border , pfcol=colors_in , plwd=1 , plty=1,
              cglcol="grey", cglty=1, axislabcol="grey", caxislabels=seq(0,1,5), cglwd=0.8,
              vlcex=0.8,
              title=paste0("Tamaño:",tamanyo))

El cuarto grupo pertenece a los habitantes propios de un aera con una gran densidad de población, sim embargo con un número de habitantes de región menor que la media, deben ser ciudades pequeñas pobladas, sus ingresos y tamaño familiar están por debajo de la media, lo llamaremos “areas urbanas densas y pequeñas”. Este el cluster con un menor número de pertenecientes.

  tamanyo<-centroidesJerarquico[5,2]
  
  radarchart( as.data.frame(centroidesOptimizacionParaRadar[c(1:3,8),])  , axistype=1 ,
              pcol=colors_border , pfcol=colors_in , plwd=1 , plty=1,
              cglcol="grey", cglty=1, axislabcol="grey", caxislabels=seq(0,1,5), cglwd=0.8,
              vlcex=0.8,
              title=paste0("Tamaño:",tamanyo))

El quinto cluster tiene todo por encima de la media, especialmente el nivel de ingresos y el tamaño familiar, dene ser ciudades con un areas grandes, con mucho trabajo y un PIB per capita muy alto, lo llamaremos “familias enriquecidas de grandes ciudades”. Cabe destacar que este cluster también cuenta con un gran número de pertenecientes.

  tamanyo<-centroidesJerarquico[6,2]
  
  radarchart( as.data.frame(centroidesOptimizacionParaRadar[c(1:3,9),])  , axistype=1 ,
              pcol=colors_border , pfcol=colors_in , plwd=1 , plty=1,
              cglcol="grey", cglty=1, axislabcol="grey", caxislabels=seq(0,1,5), cglwd=0.8,
              vlcex=0.8,
              title=paste0("Tamaño:",tamanyo))

A este grupo pertenecen habitantes de areas poco pobladas y con una baja densidad de población, el tamaño familiar es también bastante bajo mientras el nivel de ingresos esta bastante por encima de la media, por eso creo que sona areas rurales. Lo llamaremos “areas rurales enriquecidas”.

  tamanyo<-centroidesJerarquico[7,2]
  
  radarchart( as.data.frame(centroidesOptimizacionParaRadar[c(1:3,10),])  , axistype=1 ,
              pcol=colors_border , pfcol=colors_in , plwd=1 , plty=1,
              cglcol="grey", cglty=1, axislabcol="grey", caxislabels=seq(0,1,5), cglwd=0.8,
              vlcex=0.8,
              title=paste0("Tamaño:",tamanyo))

El último grupo se caracteriza por una densidad de población y número de habitantes ligeramente por debajo de la media, mientas que el nivel de ingresos y el tamaño de la unidad familiar es bastante alto, lo llamaremos “areas familiares enriquecidas”.

Por último, vamos a mapear nuestrso resultados para tener una visión global de cómo se encuentra distribuida los habitantes de cada cluster.

library(leaflet)
pal<-colorFactor(palette="Set1",
                 domain=census$cluster)
popup<-paste0('<b>tamañofam:</b> ', as.character(census$MeanHHSz), '<br>',
              '<b>ingresos:</b>', as.character(census$MedHHInc), '<br>',
              '<b>densidad:</b>', as.character(census$RegDens), '<br>' ,
              '<b>población:</b>', as.character(census$RegPop))
leaflet(data=census)%>%
  setView(lng=mean(census$LocX), lat=mean(census$LocY),zoom=2)%>%
  addTiles(urlTemplate = 'http://{s}.basemaps.cartocdn.com/dark_all/{z}/{x}/{y}.png')%>%
  addCircles(~LocX,~LocY,popup = popup,color=~pal(cluster))